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DIFFERENTIAL  EQUATIONS  WITH  APPLICATIONS 
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Grant  No.  AFOSR-75-2818 

Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

Prepared  by  Joseph  E.  Flaherty 


During  the  course  of  this  grant  algorithms  were  developed 
and  analyzed  for  the  numerical  solution  of  singularly-perturbed 
(or  stiff)  initial  value  (cf .  [4] ,  [8] ,  and  [9])  and  boundary  value 
Q-JcO^-r  [ITT  EdrdH  problems  for  ordinary  differential  equations  and 
initial-boundary  value  problems  for  partial  differential  equations 
[11] .  These  general  purpose  methods  have  been  applied  to  a  wide 
variety  of  problems  arising  in  several  disciplines,  including 


optimization  and  optimal  control  :^(ef . 


nonlinear  oscil¬ 


lations  ‘(cf  .  [Sih; — [7]-),  chemical  reactions  and  hydrodynamic 

stabilitylTeK  ‘ ^  ■  ;,  ■  ... 


Five  graduate  students  were  partially  supported  by  this' 
grant  and  assisted  Professor  Flaherty  in  his  investigations.  Four 
of  the  students  earned  MS  degrees  in  either  Applied  Mathematics  or 
Computer  Science  and  one  earned  a  Ph.D.  degree  in  Mathematics.  A 
paper  based  on  W.  Mathon's  MS  thesis  [10]  has  been  accepted  for 
publication  in  the  SIAM  J.  Sci.  and  Stat.  Comp,  and  it  will  appear 
in  print  this  Fall.  Two  other  MS  students  (G.  M.  Heitker  and  J.  M. 
Coyle)  extended  Flaherty  and  Mathon's  work  to  systems  of  boundary 
value  problems  and  boundary  value  problems  having  rapidly  oscillat¬ 
ing  solutions.  Papers  based  on  their  theses  are  being  prepared  for 
publication.  Both  Heitker  and  Coyle  are  continuing  their  studies 
towards  Ph.D.  degrees.  S.  F.  Davis  recently  completed  his  Ph.D 
dissertation  [11]  on  adaptive  grid  finite  element  methods  for  time 
dependent  partial  differential  equations.  His  method  is  capable 
of  adapting  the  finite  element  mesh  in  order  to  follow  such  sharp 


See  the  list  of  References  and  Abstracts  following  this  report. 


transitions  as  boundary  layers,  shock  layers,  wave  fronts  and  melt¬ 
ing,  boiling,  or  reaction  interfaces.  A  first  paper  based  on  part 
of  Dr.  Davis^'  dissertation  is  in  its  final  stages  of  preparation 
and  will  be  submitted  for  publication  shortly.  Additional  papers 
will  follow. 

Dr.  P.  T.  Boggs  of  the  U.  S.  Army  Research  Office  was  also 
supported  for  one  month  during  1975  and  published  a  paper  on  opti¬ 
mization  techniques  [3]  due  in  part  to  this  support. 

This  grant  has  also  provided  support  for  Professor  R.  E. 
O'Malley,  Jr.  of  the  University  of  Arizona  to  visit  R.P.I.  for  about 
one  month  each  year  and  collaborate  with  Professor  Flaherty.  This 
collaboration  has  been  most  fruitful  and  rewarding.  Five  joint 
papers  have  been  published  (cf.  [1],  [2],  [4],  [8],  [9])  and  two 
more  are  being  prepared  for  publication. 

During  the  course  of  this  investigation  Professor  Flaherty 
was  invited  to  lecture  on  material  pertaining  to  this  grant  at  the 
following  institutions:  The  Weizmann  Institute  of  Science,  Rehovot, 
Israel  (1974) ,  Clarkson  College  (1975) ,  Drexel  University  (1976) , 
the  University  of  Utah  (1976) ,  the  University  of  Melbourne,  Melbourne, 
Australia  (1977) ,  the  University  of  Arizona  (1977) ,  the  Catholic 
University,  Nijmegen,  The  Netherlands  (1978) ,  Stanford  University 
(1980),  the  University  of  British  Columbia  (1980)  and  the  University 
of  California  at  Davis  (1980) . 
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THE  NUMERICAL  SOLUTION  OF  BOUNDARY  VALUE  PROBLEMS  FOR  STIFF 

DIFFERENTIAL  EQUATIONS 

by 
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Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

and 

R.  E.  O'Malley,  Jr. 

Department  of  Mathematics 
University  of  Arizona 
Tucson,  Arizona  85721 

Abstract 

The  numerical  solution  of  boundary  value  problems  for  cer¬ 
tain  stiff  ordinary  differential  equations  is  studied.  The 
methods  developed  use  singular  perturbation  theory  to  construct 
approximate  numerical  solutions  which  are  valid  asymptotically; 
hence,  they  have  the  desirable  feature  of  becoming  more  accurate 
as  the  equations  become  stiffer.  Several  numerical  examples 
are  presented  which  demonstrate  the  effectiveness  of  these  method 

Published  in  Mathematics  of  Computation,  Vol.  31,  No.  137, 

January  1977,  pp.  66-93. 


ON  COMPUTING  SINGULAR  OPTIMAL  CONTROLS 


BY 

Joseph  E.  Flaherty 
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Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

and 

R.  E.  O'Malley,  Jr. 

Department  of  Mathematics 
University  of  Arizona 
Tucson,  Arizona  85721 

Abstract 

We  consider  singular  optimal  control  problems  consisting 
of  a  state  equation 

X  =  Ax  +  bu 

for  vectors  x  and  gcalars  u  and  a  cost  functional 
T 

if  2  2 

J  =  J  (x'Qx  +  e  u  )dt 
0 

to  be  minimized  for  |u|£  m  and  e  =  0.  By  considering  the  problem 
as  e  ->•  0 ,  singular  perturbation  concepts  can  be  used  to  compute 
solutions  consisting  of  bang-bang  controls  followed  by  singular 
arcs.  The  procedure  further  develops  a  ntimerical  technique 
proposed  by  Jacobson,  Gershwin,  and  Lele,  as  well  as  additional 
analytic  methods  of  other  authors . 

Published  in  the  IEEE  Transcations  on  Automatic  Control, 

Vol.  AC-22,  No.  4,  August  1977,  pp  640-648. 


An  Algorithm  Based  on  Singular  Perturbation  Theory, 

For  Ill-Conditioned  Minimization  Problems 

by 

Paul  T.  Boggs 

U.  S.  Army  Research  Office 
P.  0.  Box  12211 

Research  Triangle  Park,  N.  C.  27709 

Abstract 

The  problem  of  finding  an  unconstrained  local  minimum  of  a 
function  f  is  considered.  Ill-conditioned  minimization  problems 
are  characterized  by  long,  very  narrow  valleys  which  tend  to  trap 
many  numerical  procedures,  i.e.,  methods  applied  to  such  problems 
tend  to  "jam".  To  overcome  this  difficulty,  the  problem  is  re¬ 
formulated  as  that  of  integrating  the  differential  equation 
x'  =  -Vf(x).  The  ill-conditioning  of  the  minimization  problem 
carries  over  to  the  differential  equation  by  making  it  stiff. 

Thus  methods  for  integrating  stiff  systems  of  equations  are  in¬ 
vestigated  and  adaped  to  this  special  problem.  In  particular,  a 
recent  jnethod  of  Miranker  based  on  singular  perturbation  theory 
is  modified  and  shown  to  be  applicable  to  this  problem.  Some 
numerical  results  are  presented  which  indicate  that  the  resulting 
algorithms  are  quite  competitive. 

In  SIAM  J.  Numer.  Anal.,  Vol.  14,  No.  5,  Oct.  1977,  pp.  830-843. 


SINGULAR  SINGULAR-PERTURBATION  PROBLEMS 


by 

R.  E.  O'Malley,  Jr. 

Department  of  Mathematics 
^  University  of  Arizona 

Tucson,  Arizona  85721 

and 

J.  E.  Flaherty 

Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

Abstract 

We  consider  the  nonlinear  system 

(1)  e  z  =  f  (z,  t ,  e) 

of  M  nonlinear  equations  on  some  finite  subinterval  of  t  ^  0 
subject  to  the  initial  condition 

(2)  z(0)  =  z^(e) 

in  the  limit  as  the  small  positive  parameter  e  ->■  0.  Any  acquaint¬ 
ance  with  singular  perturbations  would  lead  one  to  expect  that  the 
unique  solution  of  this  initial  value  problem  might  converge  as 
e  0  to  a  solution  of  the  limiting  equation 

(3)  f(Zg,t,0)  =  0 

away  from  any  "boundary  layer”  regions  of  nonuniform  convergence. 
In  the  singular  situation  when  is  a  singular  matrix, 

this  limiting  equation  has  an  infinite  number  of  solutions,  so  the 
limiting  equation  (3)  is  not  adequate  for  the  purpose  of  determini 
the  reduced  problem  satisfied  by  the  limiting  solution  (assuming 
that  such  a  limit  exists) . 


-2- 


We  consider  the  initial  value  problem  (1)  -  (2)  under  the 
primary  assumptions  that  the  matrix 

(4)  f^(z,t,0)  fQ^(z,t) 

has  constant  rank  k,  0  <  k  <  M,  that  its  null  space  is  spanned 

by  M  -  k  linearly  independent  eigenvectors ,  and  that  f  has  k 

z 

stable  eigenvalues  (counting  multiplicities)  for  all  z  and  t. 
Assuming  sufficient  differentiability  of  f  and  asymptotic 
expansions  for  both  f  and  z^ ,  we  construct  an  asymptotic 
solution  of  the  form 

(5)  z(t,e)  =  Z{t,e)  +  n(T,£) 

where  Z  and  II  have  asymptotic  power  series  expansions  in  e 
the  terms  of  n  all  tend  to  zero  as  the  stretched  Variable 
T  =  t/e  tends  to  infinity. 


Published  in  Lecture  Notes  in  Mathematics,  No. 594, 

Singular  Perturbations  and  Boundary  Layer  Theory,  Lyon  1976, 
A.- Do Id  and  B.  Eckman ,  Eds.,  Springer -Verlag ,  Berlin,  1977, 
pp  422-436. 


FREQUENCY  ENTRAINMENT  OF  A  FORCED  VAN  DER  POL  OSCILLATOR 

by 

J.  E.  Flaherty 

Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

and 

F.  C.  Hoppensteadt 

Courant  Institute  of  Mathematical  Sciences 
New  York  University 
New  York,  New  York  10012 

ABSTRACT 

A  van  der  Pol  relaxation  oscillator  that  is  subjected  to 
external  sinusoidal  forcing  can  exhibit  stable  and  unstable 
periodic  and  almost  periodic  responses.  For  some  forcing 
amplitudes  it  even  happens  that  two  stable  subharmonics 
having  different  periods  may  coexist.  We  investigate  here  the 
stable  responses  of  such  forced  oscillators.  By  numerically 
computing  the  rotation  nxamber  of  stable  oscillations  for 
various  values  of  the  forcing  amplitude  and  oscillator  tuning, 
we  obtain  descriptions  of  regions  of  phase  locking,  successive 
bifurcation  of  stable  subharmonic  and  almost  periodic 
oscillations , and  overlap  regions  where  two  distinct  stable 
oscillations  can  coexist. 

Published  in  Studies  in  Applied  Mathematics  ,  Vol .  58,  1978 
pp  5-15. 


NTTMEJlICaL  INVESTIGATION  OF  THE 


EFFECT  OF  A  CORIOLIS  FORCE  ON  THE  STABILITY  OF  PLANE 

POISEUILLE  FLOW 

by 

J-E.  Flcihexty  and  R.  C.  D±Pi±ma 

Departmen-t  of  M^^thematical  Sciences 
Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

ABSTRACT 

.  The  staJbd-Lity  of  the  viscous  flow  between  two  parallel 
horizoirtal  plates  due  to  a  constant  reduced  pressure  gradient 
in  a  system  rotating  about  a  vertical  axis  is  studied.  The 
critical  value  of  the  Reynolds  number  R,  based  on  the  reduced 
pressure  gradient,  is  a  function  of  a  dimensionless  rotation 
parameter  T,  the  Taylor  nundjer .  A  numerical  solution  of 
the  eigenvalue  problem  shows  that  (i)  the  minimum  point 
for  the  viscous  instability  mode  associated  with  plane 
Poiseuille  flow  at  T  =  0  disappears  at  a  value  of  T  of  about 
T-0.4,  and  (ii)  for  T^O  a  new  instability  mode  appears  as 
a  result  of  the  Coriolis  effect  on  the  basic  flow  and  in 
the  perturbation  equations.  This  new  instability  gives 
the  critical  value  of  R  for  values  of  T  as  small  as  0.06. 


Published  in  The  Physics  of  Fluids,  Vol.  21,  No.  5,  May  1978, 
pp  718-726 


NONUNIQUE  STABLE  RESPONSES  OF  EXTERNALLY 
FORCED  OSCILLATORY  SYSTEMS 

by 
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Salt  Lake  City,  Utah  84112 

and 

J.E.  Flaherty 

Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Institute 
Troy,  New  York  12181 

ABSTRACT 

A  population  model,  based  on  a  system  of  difference  equations, 
is  presented  which  possesses  a  stable  relaxation  oscillation.  When 
this  system  is  subjected  to  external  periodic  forcing,  it  exhibits 
phase  locking  of  subharmonics  and  for  some  parameter  values  the 
co-existence  of  distinct  stable  periodic  solutions.  Successive 
bifurcations  of  subharmonics  can  be  described  in  terms  of  their  ro¬ 
tation  numbers. 

In  Annals  of  the  New  York  Academy  of  Sciences,  Vol.  316,  Bifurcation 
Theory  and  Applications  in  Scientific  Disciplines,  O.  Gurel  and 
O.  E.  Rossler,  Eds.,  1979. 


THE  NUMERICAL  SOLUTION  OF  SINGULAR 
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Polytechnic  Institute,  Troy,  New  York  12181,  U.S.A. 

R.E.  O'Malley,  Jr. 

Department  of  Mathematics  and  Program  in  Applied 
Mathematics,  University  of  Arizona,  Tucson, 
Arizona  85721,  U.S.A. 


ABSTRACT 


We  consider  the  vector  initial  value  problem  ey  =  f(y,t,e), 
y(0)  -  y^(e)  in  the  situation  when  the  m  x  m  matrix  fy(y,t,0)  is 
singular  with  constant  rank  k  <  m  and  has  k  stable  eigenvalues. 

We  show  how  to  determine  the  unique  limiting  solution  Y^  of  the 
reduced  problem  f(YQ,t,0)  =  0  and  how  to  obtain  a  uniform  asymptotic 
expansion  of  the  solution  which  is  valid  for  small  values  of  e 
on  finite  t  intervals.  A  numerical  technique  is  developed  to 
calculate  the  limiting  solution  and  the  results  of  some  examples 
are  compared  with  an  existing  code  for  stiff  differential  equations. 


In  Numerical  Analysis  of  Singular  Perturbation  Problems,  P.  W. 
Hemker  and  J.  J.  H.  Miller,  Eds.,  Academic  Press,  London,  1979, 
pp.  143-180. 
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SINGULAR  SINGULARLY-PERTURBED  INITIAL  VALUE  PROBLEMS 
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and 

J.  E.  Flaherty 
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Troy,  New  York  12181 

ABSTRACT 

This  paper  considers  the  vector  problem  ey  =  f(y,t,e),  y(0) 

prescribed,  when  the  m  m  matrix  fy(y,t,0)  is  singular  and  of 

rank  k,  0  ^  k  <  m,  with  k  stable  eigenvalues.  The  reduced  system 

f(YQ,t,0)  =  0  is,  then,  insufficient  to  determine  the  limiting 

solution  Yq  for  t  >  0.  Instead,  depending  on  the  row-reduced 

form  of " f  ,  a  differential  equation  is  obtained  for  Y„  and  an 

y  0 

initial  vector  ^^(0)  is  obtained  as  the  steady-state  solution  of 
an  inner  problem.  Alternative  analytical  and  numerical  possi¬ 
bilities  are  explored,  and  analytical  and  numerical  solutions 
are  obtained  for  representative  problems. 

In  SIAM  J .  Appl.  Math. ,  Vol.  38,  No.  2,  April  1980,  pp.  225-248. 


^  COLLOCATION  WITH  POLYNOMIAL  AND  TENSION 
SPLINES  FOR  SINGULARLY  PERTURBED 


BOUNDARY  VALUE  PROBLEMS* 


Joseph  E.  Flaherty 
Department  of  Mathematical  Sciences 
Rensselaer  Polytechnic  Institute 

and 

William  Mathon 
IBM 

Federal  Systems  Division 


ABSTRACT 

Collocation  methods  using  both  cubic  polynomials  and  splines 
in  tension  are  developed  for  second  order  linear  singularly- 
perturbed  two-point  boundary  value  problems.  Rules  are  developed 
for  selecting  tension  pareimeters  and  collocation  points.  The 
methods  converge  outside  of  boundary  layer  regions  without  the 
necessity  of  using  a  fine  discretization.  Numerical  examples 
comparing  the  methods  are  present. 
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ABSTRACT 


A  finite  element  method  is  developed  to  solve  initial¬ 
boundary  value  problems  for  vector  systems  of  parabolic  partial 
differential  equations  in  one  space  dimension  and  time.  The 
methods  adapt  the  computational  mesh  as  the  solution  progresses 
in  time  and  are  thus  able  to  follow  and  resolve  such  sharp  tran¬ 
sitions  as  boundary  layers,  shock  layers,  or  wave  fronts.  This 
permits  an  accurate  solution  to  be  calculated  with  fewer  mesh 
points  than  would  be  necessary  with  a  uniform  mesh. 

The  overall  method  contains  two  parts,  a  solution  algorithm 
and  a  mesh  selection  algorithm.  The  solution  algorithm  is  a 
finite  element-Galerkin  method  on  trapezoidal  space- time  elements, 

using  either  piecewise  linear  or  cubic  polynomial  approximations. 

(continued) 


An  efficient  and  stable  algorithm  is  developed  to  solve  the  re¬ 
sulting  block  tridiagonal  nonlinear  discrete  system. 

The  mesh  selection  algorithm  builds  upon  the  work  of  de 
Boor,  Jupp  and  others  on  variable  knot  spline  interpolation. 

Both  the  solution  and  the  mesh  selection  algorithm  are 
analyzed  in  detail .  Error  estimates  are  obtained  for  the  solution 
algorithm.  The  mesh  selection  algorithm  is  shown  to  produce  a 
mesh  which  minimizes  the  L2  error  of  the  solution  for  a  given 
nxamber  of  mesh  points. 

A  computer  code  implementing  these  algorithms  has  been 
written  and  applied  to  a  number  of  problems.  These  computations 
confirm  that  the  theoretical  error  estimates  are  attained  and 
demonstrates  the  utility  of  variable  mesh  methods  for  partial 
differential  equations. 


